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Abstract 

We present a brief analysis of the crossover phase diagram for the decay of a 
metastable phase in a simple dynamic lattice-gas model of a two-phase system. We 
illustrate the nucleation-theoretical analysis with dynamic Monte Carlo simulations of 
a kinetic Ising lattice gas on square and cubic lattices. We predict several regimes 
in which the metastable lifetime has different functional forms, and provide esti- 
mates for the crossovers between the different regimes. In the multidroplet regime, 
the Kolmogorov-Johnson-Mehl-Avrami theory for the time dependence of the order- 
parameter decay and the two-point density correlation function allows extraction of 
both the order parameter in the metastable phase and the interfacial velocity from the 
simulation data. 
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1 Introduction 



The decay of a metastable phase through nucleation and growth of droplets of a more stable 
phase is a common feature of phase-transformation processes in both natural and techno- 
logical settings. It occurs both in systems where the phases involved are non-crystalline 
and crystalline Often such processes are well described by a combination of nucleation 
theory [|2|] and the Kolmogorov-Johnson-Mehl-Avrami (KJMA) theory [J3J, U 0], which de- 
scribes nucleation followed by irreversible growth of the nucleated droplets |@]. 

Here we present some recent theoretical results on the effects of finite system size in 
simple model systems with microscopic dynamics, which undergo decay of a metastable 
phase through homogeneous nucleation and growth. These results are illustrated and con- 
firmed by extensive dynamic Monte Carlo (MC) simulations. The simulations utilize sev- 
eral advanced algorithms [[7|, ^ [TT], [11]] that enable numerical study of extremely long- 
lived metastable phases without changing the underlying dynamics. Although the mod- 
els that we study are highly simplified, there is evidence from magnetic systems that the 
predicted finite-size effects are experimentally observable [[I2|, [T3| ]. Moreover, our data- 
analysis methods should also be applicable to experiments. 

The rest of this paper is organized as follows. The theoretical model is presented in 
Sec. 0, and relevant aspects of the nucleation and KJMA theories are presented in Sec. ^ 
Results of MC simulations are given in Sec. 0. In particular, we present a "metastability 
phase diagram," metastable lifetimes for two- and three-dimensional systems, and measure- 
ments of the metastable order parameter and the propagation velocities of domain walls in 
two-dimensional systems. A discussion and conclusions are given in Sec. §. 



2 Theoretical Model 

We consider a simple lattice-gas model of a two-phase system, defined by the effective 
Hamiltonian 

TihG = -<P C * C j ~ A* H °i ■ CD 

Here qg{0, 1} is a local concentration variable at the ith site of a <i-dimensional hypercubic 
lattice with periodic boundary conditions, (f> is an attractive interaction between neighboring 
sites, and \x is a chemical potential favoring q=1. The sums J2(i,j) and J2i mn over ai l 
nearest-neighbor pairs and all lattice sites, respectively, and the lattice constant is taken as 
unity. The order parameter conjugate to /x is the dimensionless density (for d = 2 often 
called "coverage"), 9 = A/" -1 Y,i c u where J\f is the total number of sites. Below a critical 
temperature T c , two uniform phases coexist along a line of first-order phase transitions at 
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p = po = —dcj). Here p can be interpreted as an (electro)chemical potential, or related to 
an (osmotic) pressure p as p — /x = ln(p/p ) where fc B is Boltzmann's constant, T the 
absolute temperature, and p the pressure at coexistence. Analogous relations can be made 
to supersaturation, supercooling, etc. Since the distance of p from the coexistence curve is 
the important parameter, we introduce the notation Ap = p — p Jl4|]. 

The phase degeneracy that leads to coexistence is lifted for Ap ^ 0. For Ap > the 
high-density phase with 9 near unity is the equilibrium phase, while Ap < favors the 
low-density phase. In both cases the uniform phase which is not favored by the sign of Ap 
can be extremely long-lived. The kinetics of the decay towards equilibrium of a system 
initially in such a long-lived metastable phase is the topic of this paper. 

The lattice-gas model does not have an intrinsic dynamic. However, on time scales 
that are long compared to an inverse phonon frequency, the approach to equilibrium can 
be described by a stochastic dynamic. Here we study local dynamics, in which the density 
is not conserved. Among the many possible dynamics that satisfy detailed balance, and 
thus ensure approach to equilibrium, we choose the single-site Glauber dynamic [[15]] with 
updates at random sites. In this dynamic a proposed update of a single q is accepted 
with the probability W(ci— >1— Ci) = [l+exp (AE / ksT)}' 1 , where AE is the energy change 
which would ensue if the update were accepted. Each attempted update advances the clock 
by one step, and we measure time in MC Steps per Site (MCSS). 

3 Nucleation and Growth 

Next we summarize elements of the theory of phase transformation by nucleation and 
growth in systems with nonconserved order parameter. For further details, see Refs. [|J 
[l2t |16| ]. We concentrate on homogeneous nucleation, thus ignoring effects of surfaces or 
disorder. Hence our comparisons will be with simulations using periodic boundary condi- 
tions and position-independent interactions and chemical potentials. 

A central role is played by the metastable lifetime, (r(Ap, T, L)), which is the typical 
time it takes to convert half of the system volume from the metastable to the stable phase. 
For concreteness we define it as the mean first-passage time to 9 = 1/2 from an initial 
condition with all q =1 and Ap < 0. 

For the class of attractive lattice-gas models with short-range interactions studied here, 
the fluctuations that lead to the decay of the metastable phase are compact, (/-dimensional 
droplets of radius R. The free energy of the droplet has two competing terms: a positive 
interface term oc and a negative bulk term oc \Ap\R d . The competition between 
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these yields a critical droplet radius, 

R c (A(i,T) 



\Afi\Ae (T) 



(2) 



where er(T) is the equilibrium surface tension and A9 (T) is the equilibrium density dif- 
ference at coexistence. Both are positive for T < T c . Droplets with R < R c tend to decay, 
while droplets with R > R c tend to continue growing. The free-energy cost of the critical 
droplet with R = R c , relative to the uniform metastable phase, is 



AF(Afi,T) = n d a(T)' 



d-1 



(|A/i|A0oCr) 



d-1 



(3) 



where Q d is a weakly T-dependent shape factor such that the volume of a droplet of radius 
R is VtdR d - Nucleation is a stochastic process, and the nucleation rate per unit volume is 
given by a Van't Hoff-Arrhenius relation, 



J(A/i,T) oc |A/i| K exp 



S(T) 



k B T\A^ 



(4) 



where S(T) is the A/^-independent part of AF. For the model studied here, the prefactor 
exponent K equals 3 for d = 2 [§, [R], [17|, 0] and is believed to be -1/3 for d = 3 10]. 

For systems of finite linear size L, a T- and L-dependent crossover called the Thermo- 
dynamic Spinodal (THSP), | A/i| TH sp [ |16H is determined by the condition R c w L. Thus 
IA^Ithsp ~ L^ 1 . For |A//| < | A/i| TH sp, R c would exceed L. This is called the Coexis- 
tence (CE) regime because the critical fluctuation here resembles two coexisting domains 
of different density. For \Afi\ > |A/i| T hsp (but not too large, as discussed below), the 
lifetime is dominated by the inverse of the total nucleation rate: 



r SD (A//,T,L)) « (L /(A/x, T) ) oc L' d \ Av\~ K exp 



S(T) 



^TIA^I^- 1 



(5) 



This is inversely proportional to the system volume, L d . The subscript SD stands for Single 
Droplet and indicates that the phase transformation is completed by the first critical droplet 
which nucleates. 

The fact that a supercritical droplet grows at a finite velocity, t>(A/i,T), leads to a 
second crossover, called the Dynamic Spinodal (DSP) [|16|]. A reasonable criterion to locate 
the DSP is that the average time between nucleation events, (tsd), should equal the time it 
takes a droplet to grow to a size comparable to L. This gives the asymptotic relation 



\Afi\ mP {T,L) 



H(T) 



(d + l)k B T\nL 



(6) 
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The DSP is shown in Fig. [I] for two different values of L. The convergence to zero as 
L— >oc is exceedingly slow, so that the crossover is observable in sufficiently small sys- 
tems [[T2|, [T3p. For |A/i| > | A/i| DS p, the metastable phase decays through many droplets 
which nucleate and grow independently in different parts of the system. This is called the 
Multidroplet (MD) regime [|16|]. The KJMA theory of metastable decay in large systems 
® B Si gi yes the lifetime in this regime, 



<t md (A//,T)> 



d+1 



oc 



(d+1) In 2 



-d -K 

v d +^ \Afi\ d + 1 exp 



S(T) 



(d + l)k B T\Afi\ 



d-l 



(7) 
(8) 



which is independent of L. Note that the velocity depends on T and Afx. Comparing 
Eqs. (Q) and (||), it becomes clear that it is useful to plot hi(r) vs 1/| A/x| d_1 . In such plots, 
which are shown in Fig. 0, the lifetimes in both regimes lie on approximately straight lines, 
with the slope in the SD regime (d+1) times that in the MD regime. 

In the MD regime the order parameter decays according to "Avrami's law," 



(9(t))-9 s (Afi,T) 
9 ms (Afx,T) - 9 s (Afi,T) 



cxp 



-In 2 



(tmd) 



(9) 



where 9 s (Afi, T) and 6> ms (A/i, T) are the equilibrium and quasi-equilibrium values of the 
density in the stable and metastable phases, respectively. Equilibrium simulations or mea- 
surements easily yield # s (AyU, T). By fitting Eq. @ to simulational or experimental data 
for 9(t), one can therefore estimate # ms (A/i, T) and (t M d(A^, T)). The latter gives the 
combination Iv d of the nucleation rate and growth velocity [Eq. (J7J)], but not / and v sepa- 
rately. These quantities can, however, be resolved by using results for the two-point density 
correlation functions in the KJMA picture [ |T9| , |2D| ] . These give the time-dependent variance 



(9' 2 (t)) - (9(t)) 2 « L' d dn d {(9(t)) - 9 S } Z §(Iv d t d+1 ) (2t) d v 



(10) 



where $ is determined by numerical integration of an expression involving the correlation 
function [jl/J, [21]]. The last factor of v d in Eq. ( JTOj ) is not multiplied by /. This permits the 



separate determination of v by measurement of the variance of 9(t) in a series of indepen- 
dent simulations or experiments. 

The description of the decay in terms of droplets breaks down for sufficiently large 
\Aji\. The crossover marking the transition to this Strongly Forced (SF) regime is estimated 
as the L-independent "Mean-Field Spinodal" (MFSP), | A^| M fsp(T), for which R c = 1 /2. 
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4 Numerical Results 



The "metastability phase diagram" for the square-lattice Ising lattice-gas model is shown 
in Fig. III. The three decay regimes shown are the SF, MD, and SD regimes. The CE regime, 
which occurs at very low |A/i|, is not shown. The MFSP, which separates the SF regime 
from the other regimes, is given by Eq. (0) with exactly known expressions for both a(T) 
and A9 . This curve goes to zero at T c . The DSP is shown for two system sizes, L=24 
and 240. Our estimate for the location of the DSP is where the standard deviation of the 
lifetime, At, is equal to (t}/2. In the SD (and CE) regimes At~(t), while Ar <C (t) in 



the MD and SF regimes [16]. 



As L increases, the portion of the diagram occupied by the MD regime increases. For 



2 < 2| A/x|/0 < 4, the critical droplet consists of a single lattice site Q22| l. For the dynamic 
used here, there exists an analytical low-temperature estimate for the DSP in this range of 
| A/i| [[23]] . It is obtained by equating the growth time and the nucleation time, and is (for 
d = 2) given by 

2|A,| DSP (r,L) =4 _4j^/3 lnL _ N . (ID 

<p V2 / 

The only adjustable parameter is the last term inside the parenthesis, which was fit to give 
agreement for L between 8 and 240. This expression for the DSP is marked by the two 
solid straight lines in Fig. [[]. It separates the SF regime from a low-temperature SD regime, 
where the nucleating droplet occupies only a single lattice site. The agreement between 
the simulation data and the analytical expression is very good. Below the intersection with 
the MFSP, the low-temperature form for | A/i|dsp merges smoothly into the form given by 
Eq. (g), which goes to zero as T — ► T c . 

As seen fromEqs. (|5|) and (|8|), plotting (r) versus 1/| A/z| d_1 gives a reasonable estimate 
for S(T) /k B T in the SD regime and H(T) /(d+ l)k B T in the MD regime. Examples are 
shown in Fig. ^|for d=2 (a) and d=3 (b). Slight curvatures are present, due to prefactors in 
the nucleation rate and also to the proximity of a crossover (the THSP, DSP, or MFSP). Well 
away from crossovers, the curvatures have been used to estimate the prefactor exponent 
K [|5J, [R]]. The symbols in both panels (with error estimates from Ar), are from 1000 
simulations using the n-fold way algorithm [J7|]. 

The ra-fold way algorithm is based on dividing all sites into classes reflecting the current 
state of the site and its neighbors. It yields class populations of n=10 classes for d=2 and 
n=14 classes for d=3. These populations can be utilized to 'coarse-grain' or 'lump' all 
configurations onto the one-dimensional variable 9. Measuring these class populations 
and utilizing the analytical form for the acceptance probabilities yields average growth and 
shrinkage rates for the stable phase. This is the central idea behind the Projective Dynamics 
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(PD) method [§, |9|, [H]], which was used to obtain the curves in Fig. f2[ Here we measured the 
class populations in equilibrium, which provides reasonable approximations for the actual 
class populations during metastable escape for small |A/i|. This approximation causes 
the slight differences between the curves and the measured data points. Using the correct 
measured class populations during the escape from the metastable phase gives the exact 
lifetime (within statistical errors) [[TOtl ■ The PD method also enables one to obtain results 
for large systems from class populations in small systems [f| §]. In fact, the PD results 
for the large systems in Fig. [| are obtained from the equilibrium class populations in the 
smaller systems. By using such advanced algorithms [|7|, D, [LOj, [TT|], one can speed up 
calculations of metastable lifetimes by many orders of magnitude — without changing the 
dynamics of the system. 

The approximate KJMA theory of phase transformation kinetics is a remarkably ver- 
satile tool for analysis of experimental results. Nevertheless, there are very few stud- 
ies in which it has been explicitly tested in theoretical models with microscopic dynam- 
ics on a length scale much smaller than R c ffij [24] , [Hp . Here we present results us- 
ing Eqs. @ and ( |10| ) to analyse the decay of the metastable phase in d=2 MC simula- 
tions in the MD regime at 0.8T C [[HJ|. The systems are sufficiently large to contain at 
least several hundred supercritical droplets (L=256 for 2|A/i|/0 > 0.20 and L=1024 for 
0.20 > 2\Afx\/(f) > 0.12). 

Figure |3] shows the metastable density 9 ms (Afi, T) for A/i < 0, obtained by fitting 
data for (9(t)) to Eq. (Q). The equilibrium density # s (A/a, T) is shown for A/i > 0. As 
one would expect, the metastable phase becomes gradually less well defined as |A/i| is 
increased. This is reflected by the increasing differences between the data points shown as 
empty circles and squares, which were obtained using different data cutoffs in the fitting 
procedure [[21)] . For smaller | A/i| these estimates coincide, indicating that the details of the 
fitting procedure become unimportant as the metastable phase becomes more well defined. 
As a check on the procedure, metastable densities measured by a transfer-matrix method 
[ ]2~T1 ] are also shown. The agreement between the estimates obtained by these two unrelated 
methods is remarkable over the whole range of | A/i| shown. 

In Fig. [| we show the average velocity of the convoluted interface between the metastable 
and stable phases, v (A/i, T), obtained by fitting Eq. ( JT0| ) to the variance of the simulated 
densities, while using the values of 9 ms and Iv 2 obtained from fitting (9(t)) to Eq. (Q). 
As a check on this rather indirect way of estimating the velocity, the figure also contains 
data from a simulation of a flat interface propagating into a metastable bulk phase in which 
nucleation is suppressed [|2~6|], and a curve obtained from an analytical nonlinear-response 
approximation [[TTJ, . The agreement between the different estimates is good over the 
whole range of | Afx\ shown. 
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5 Discussion and Conclusion 

We have shown that even for the simple case of homogeneous nucleation and growth appli- 
cable to the decay of the metastable phase in a nearest-neighbor Ising lattice-gas model with 
periodic boundary conditions, there are different decay regimes separated by crossovers 
(Fig. [I]). This is because there are four length scales in the problem. Two length scales 
are fixed, the lattice spacing and the linear system size L. Two length scales depend on T 
and the distance |A//| from the coexistence curve. These are the critical droplet radius R c 
and the typical size to which a supercritical droplet grows before it interacts with another, 
Ro ~ (tmd)v- Starting from the smallest |A/x|, the different decay regimes are the CE 
regime where 2R C > L, the SD regime where 2R c <t^L<^R and a single nucleating droplet 
is responsible for the decay of the metastable phase, the MD regime where Rq<^L and 
many supercritical droplets are formed during the decay process, and the SF regime where 
the droplet picture is not applicable. These regimes are separated by crossover lines, which 
depend on T and L in different ways. Between the CE and SD regimes is the Thermody- 
namic Spinodal with |A/i|THSP ~ 1/L. Between the SD and MD regimes is the Dynamic 
Spinodal with |A/x|dsp ~ l/(lnL)~ asymptotically for very large L. Between the MD 
and SF regimes is the Mean-Field Spinodal, which is independent of L. 

The average metastable lifetime (r) takes different functional forms in the different 
regimes (Fig. 0). For example, in the SD regime it is inversely proportional to the volume 
[Eq. (]5J)], while in the MD regime it is independent of the volume [Eq. (Q)]. In the MD 
regime, "Avrami's law" for the time dependence of the order-parameter decay yields 9 ms 
and the combination Iv d of the nucleation rate / and the interface velocity v. However, 
including two-point density correlation functions enables one to obtain / and v separately. 
By fitting the KJMA theory to simulation data for d=2 we obtain results for the metastable 
density (Fig. |]) and interfacial velocity (Fig. 0) that agree with unrelated methods of cal- 
culating the same quantities. We believe the data-analysis methods discussed here should 
be applicable to many experimental systems, as well. 
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Figure 1: "Metastability phase diagram" for d=2, with the Single Droplet (SD), Multidroplet 
(MD), and Strongly Forced (SF) regimes labeled. The coexistence (CE) regime for very small 
| A/j, | is not shown. The dot-dashed line is the estimate for the Mean-Field Spinodal (MFSP), 
as described in the text. The estimates for the Dynamic Spinodal (DSP) are shown for L=24 
(empty circles) and L=240 (filled circles). These data points are joined by dashed lines as guides 
to the eye, and agree with the analytical low-temperature expression for 2| A/z|dsp/0 > 2 from 
Eq. (11), marked by the solid straight lines. 
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Figure 2: The average lifetime of the metastable phase, measured in MCSS, shown vs 
(</>/2| A//|) a! ~ 1 for different system sizes and dimensions. In both parts, the curve for the smaller 
system is solid and the curve for the larger system is dotted. For d=2 (a), the system sizes shown 
are 20 2 and 160 2 at T=0. 3250=0. 57T C . For d=3 (b), the system sizes shown are 10 3 and 40 3 at 
T=O.50=O.44T C . The symbols are from actual simulations. The curves are PD extrapolations 
based on equilibrium data sampled for the smaller system size. See the text for further detail. 
The different slopes in the MD and SD regimes are clearly seen. After Ref. with additional 
data. 
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Figure 3: Metastable and stable densities for d=2 at 0.8T C . The equilibrium density, # s (A/x, T), 
is shown for A/x > 0. The data points represented by empty diamonds were obtained from 
equilibrium MC simulations with L=256, while the thin solid curve is from a transfer-matrix 
calculation with strip width xV=9. The metastable density, # ms (A/x, T) is shown for A/x < 0. 
The data points represented by empty circles and squares were obtained by fitting MC data 
for (6(t)) to Eq. (Q) as described in the text, while the filled circles represent transfer-matrix 
calculations with N between 5 and 9. The thin solid curve connecting the latter points is a 
guide to the eye only Statistical uncertainties in the MC results are everywhere smaller than the 
symbol size. After Ref. [21]. 
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Figure 4: The average propagation velocity of the interface separating the regions of metastable 
and stable phase for d=2 at 0.8T C . The empty squares and circles represent estimates based on 
Eqs. (Q) and ( flO| ) as described in the text. The error bars indicate the statistical uncertainties. 
The filled diamonds were obtained from MC simulations of a flat interface of length 1000, prop- 
agating into a region of metastable phase in which nucleation was suppressed [p6|]. The solid 
curve is an analytic nonlinear-response approximation [21, 2(J. After Ref. [^l|] with additional 
data. 
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